Breve introducción a la cartografía con R

IIIRqueR. Sevilla 2024

Author
Affiliation

Dr. Dominic Royé

Fundación para la Investigación del Clima (FIC)

Representación cartográfica con ggplot2

Quien Está familarizado con las posibilidades y funcionalidades del paquete ggplot2, le será fácil usar geometría espacial. Aquí ampliaremos los aspectos relacionados con la representación cartográfica utilizando algunas de esas funciones y otras nuevas. La combinación de la gramática de ggplot2 con la gestión de la información espacial que hace sf, produce un lenguaje de visualización muy intuitivo y escalable. Esta interacción se lleva a cabo principalmente con la función geom_sf(), que es aplicable a todos los tipos de geometrías en el paquete de ggplot2.

sf es compatible con funciones muy conocidas del paquete dplyr (tidyverse).

https://r-spatial.github.io/sf/reference/tidyverse.html

Para ilustrarlo con un ejemplo sencillo, cargaremos los paquetes tidyverse y sf y leemos un dataset que contiene información sobre los husos horarios. Además, cargamos otros paquetes que usaremos en esta sección.

library(tidyverse)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4          ✔ readr     2.1.5     
✔ forcats   1.0.0          ✔ stringr   1.5.1     
✔ ggplot2   3.5.1.9000     ✔ tibble    3.2.1     
✔ lubridate 1.9.3          ✔ tidyr     1.3.1     
✔ purrr     1.0.2          
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(sf)
Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE
library(RColorBrewer)
library(giscoR)
library(mapview)
library(janitor)

Adjuntando el paquete: 'janitor'

The following objects are masked from 'package:stats':

    chisq.test, fisher.test
library(mapSpain)
library(colorspace)
library(osmdata)
Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright
library(rmapshaper)
library(giscoR)
time_zone <- st_read("./data/ne_10m_time_zones.shp")
Reading layer `ne_10m_time_zones' from data source 
  `C:\Users\xeo19\Dropbox\RqueR\data\ne_10m_time_zones.shp' using driver `ESRI Shapefile'
Simple feature collection with 120 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -180 ymin: -90 xmax: 180 ymax: 90.00021
Geodetic CRS:  WGS 84

Existen múltiples funciones de operaciones para datos vectoriales.

st_crs(time_zone) # sistema de coordinadas/proyeccion
st_coordinates(time_zone) # extraer coordinadas
st_geometry(time_zone) # objeto geometry set
st_cast(time_zone, "MULTILINESTRING") # cambiar geometría
st_area(time_zone) # error
st_make_valid(time_zone) |> st_area() # calcular area

También podemos crear objetos sf con st_as_sf().

datos <- data.frame(lon = c(-3.703889, -8.533333, -5.9925, 2.166667, -3.8,
                             -1.644167, -0.9),
                       lat = c(40.4125, 42.883333, 37.3925, 41.4, 43.466667,
                             42.818333, 41.65),
                       nombre = c("Madrid", "Santiago", "Sevilla", "Barcelona",
                                "Santander", "Pamplona", "Zaragoza"))

points_sf <- st_as_sf(datos, 
                      coords = c("lon", "lat"), 
                      crs = 4326)
mapview(points_sf)

Cuando visualizamos la geometría en formato sf solo tenemos que indicarlo con la función geom_sf().

ggplot(time_zone) +
         geom_sf()

# ggplot expande el plot habitualmente
ggplot(time_zone) +
         geom_sf() +
           coord_sf(expand = FALSE)

# reproyectar on fly
ggplot(time_zone) +
         geom_sf() +
           coord_sf(crs = 3035) # ETRS89-extended / LAEA Europe

ggplot(time_zone) +
         geom_sf() +
           coord_sf(crs = "ESRI:54024") # Bonne

ggplot(time_zone) +
         geom_sf() +
           coord_sf(crs = "ESRI:54009") # Mollweide

Ejemplo: densidad de puntos

Vistas las ventajas de geom_sf(), usaremos esta función en los siguientes ejemplos. Tenemos a nuestra disposicón el paquete osmdata que nos permite crear y hacer las consultas directamente desde el entorno de R a Open Street Maps. Aún así, el uso de la overpass-turbo.eu puede ser útil cuando no estamos seguros de lo que buscamos o tenemos alguna dificultad en construir la consulta.

Antes de crear una consulta, debemos conocer qué podemos filtrar. La función available_features( ) nos devuelve un listado amplio de las características disponibles en OSM que a su vez tienen diferentes categorías (tags). Están disponibles más detalles en la wiki de OSM aquí. Por ejemplo, la característica shop contiene como categoría entre otros supermarket, fishing, books, etc.

head(available_features())
[1] "4wd_only"  "abandoned" "abutters"  "access"    "addr"      "addr:city"
head(available_tags("shop"))
# A tibble: 6 × 2
  Key   Value                                                   
  <chr> <chr>                                                   
1 shop  [[ Too many Data Items entities accessed. |  hunting  ]]
2 shop  agrarian                                                
3 shop  alcohol                                                 
4 shop  anime                                                   
5 shop  antiques                                                
6 shop  appliance                                               

En la primera parte de la consulta debemos indicar el lugar donde queremos extraer la información. La función getbb( ) crea un rectángulo de selección para un lugar dado, buscando el nombre. La función principal es opq( ) que construye la consulta final. Añadimos con la función add_osm_feature( ) nuestros criterios de filtro. Existen varios formatos para obtener el resultado de la consulta. La función osmdata_*( ) envía la consulta al servidor y en función del sufijo * sf/sp/xml nos devuelve el formato simple feature, spatial o XML.

Como ejemplo de puntos de interés (PDIs) de la base de datos de Open Street Maps (OSM), extraeremos las gasolineras.

#rectángulo de selección para la Península Ibérica
m <- c(-10, 30, 5, 46)

#construcción de la consulta
q <- m |> 
      opq(timeout = 60*100) |>
         add_osm_feature("amenity", "fuel")

str(q)
List of 5
 $ bbox     : chr "30,-10,46,5"
 $ prefix   : chr "[out:xml][timeout:6000];\n(\n"
 $ suffix   : chr ");\n(._;>;);\nout body;"
 $ features : chr "[\"amenity\"=\"fuel\"]"
 $ osm_types: chr [1:3] "node" "way" "relation"
 - attr(*, "class")= chr [1:2] "list" "overpass_query"
 - attr(*, "nodes_only")= logi FALSE
#consulta 
gas_station <- osmdata_sf(q)
gas_station
Object of class 'osmdata' with:
                 $bbox : 30,-10,46,5
        $overpass_call : The call submitted to the overpass API
                 $meta : metadata including timestamp and version numbers
           $osm_points : 'sf' Simple Features Collection with 70914 points
            $osm_lines : 'sf' Simple Features Collection with 101 linestrings
         $osm_polygons : 'sf' Simple Features Collection with 7894 polygons
       $osm_multilines : NULL
    $osm_multipolygons : 'sf' Simple Features Collection with 58 multipolygons
# limites
pi <- gisco_get_nuts(country = c("Spain", "Portugal"), nuts_level = "3", resolution = "20") |> 
          st_crop(xmin = -10, xmax = 5, ymin = 35, ymax = 45)
Warning: attribute variables are assumed to be spatially constant throughout
all geometries
bound <- group_by(pi, CNTR_CODE) |> summarise() |>
           ms_innerlines() 
plot(bound) # frontera Portugal, España

gs <- st_filter(gas_station$osm_points, pi) # filtramos las estaciones de la PI

#mapa final del resultado
ggplot(gs)+
  geom_sf(data = pi, fill = "grey90", color = "white") +
  geom_sf(colour = "#08519c",
          alpha = .1,
          size = .5) +
  geom_sf(data = bound, linewidth = .8, color = "white") +
  theme_void()

Ejemplo: Red hidrográfica

Para un ejemplo de un mapa de líneas, representaremos las redes de ríos en la cuenca del Ebro. La siguiente parte es opcional dado que hemos preparado una subset de España a partir de la versión Europea de HydroRivers (V10). Descarga de la base datos original: https://www.hydrosheds.org/products/hydrorivers

# capas disponibles en una gdb
st_layers("./data/HydroRIVERS_v10_eu.gdb") # geodatabase 

# importar la 
rivers <- st_read("./data/HydroRIVERS_v10_eu.gdb")

# límites España
esp <- esp_get_country(epsg = "4326")

# filtramos a España
rivers_esp <- st_crop(rivers, esp)

La la idea es mapear el grosor de línea y la tranparencia según el caudal medio.

# importamos los ríos de España
rivers_esp <- st_read("./data/rivers_espHydrov10.geojson") 
Reading layer `rivers_espHydrov10' from data source 
  `C:\Users\xeo19\Dropbox\RqueR\data\rivers_espHydrov10.geojson' 
  using driver `GeoJSON'
Simple feature collection with 36807 features and 15 fields
Geometry type: MULTILINESTRING
Dimension:     XY
Bounding box:  xmin: -9.472917 ymin: 36.02292 xmax: 4.320511 ymax: 44.08104
Geodetic CRS:  WGS 84
# límites de la cuenca de Ebro
hyd_ebro <- st_read("./data/cuencas_ebro.geojson") 
Reading layer `cuencas_ebro' from data source 
  `C:\Users\xeo19\Dropbox\RqueR\data\cuencas_ebro.geojson' using driver `GeoJSON'
Simple feature collection with 1 feature and 36 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -4.404174 ymin: 40.33515 xmax: 2.015186 ymax: 43.18007
Geodetic CRS:  WGS 84
# filtramos la cuenca 
rivers_esp <- st_filter(rivers_esp, hyd_ebro) 

# mapeamos el caudal medio por tamaño de línea y transparencia
ggplot(rivers_esp, 
       aes(linewidth = DIS_AV_CMS, # caudal medio
           alpha = DIS_AV_CMS)) +
  geom_sf(colour = "#2171b5") +
  scale_linewidth(range = c(.01, 3)) + # rango de grosor
  scale_alpha(range = c(.3, 1)) + # rango de transparencia
  labs(linewidth = NULL, alpha = NULL, title = expression("Caudal medio "(m^2/s))) +
  theme_void()

# expression("Caudal medio "(m^2/s)~"otro texto") ~ ó *

ggplot(rivers_esp, 
       aes(linewidth = DIS_AV_CMS, # caudal medio
           alpha = DIS_AV_CMS)) +
  geom_sf(colour = "#2171b5") +
  scale_linewidth(range = c(.01, 3), breaks = c(50, 150, 300, 500, 600)) +
  scale_alpha(range = c(.3, 1), breaks = c(50, 150, 300, 500, 600)) +
  labs(linewidth = NULL, alpha = NULL, title = expression("Caudal medio "(m^2/s))) +
  theme_void()

Ejemplo: símbolos proporcionales

Un mapa de símbolos proporcionales usa el tamañó de un símbolo, habitualmente círculos para representar una cantidad. Volvemos a usar el paro registrado por la EPA en el último semestre del año pasado a nivel provincial.

# importar datos de paro EPA INE
paro <- read_csv2("data/3996.csv") |> clean_names()
ℹ Using "','" as decimal and "'.'" as grouping mark. Use `read_delim()` for more control.
Rows: 5724 Columns: 5
── Column specification ────────────────────────────────────────────────────────
Delimiter: ";"
chr (4): Sexo, Provincias, Tasas, Periodo
dbl (1): Total

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# filtramos
paro_prov <- filter(paro, sexo == "Ambos sexos", provincias != "Total Nacional",
                    tasas == "Tasa de paro de la población", periodo == "2023T4") |> 
                dplyr::select(provincias, total) |> 
                mutate(prov_cod = str_extract(provincias, "[0-9]{2}"))

# obtenemos las provincias
prov_esp <- esp_get_prov() 

# unimos
paro_prov_sf <- left_join(prov_esp, paro_prov, by = join_by(cpro == prov_cod)) |> st_centroid()
Warning: st_centroid assumes attributes are constant over geometries
ggplot() +
  geom_sf(data = prov_esp, fill = "grey95", colour = "white") +
  geom_sf(data = paro_prov_sf, aes(size = total)) +
  labs(size = "Tasa", title = "Paro en el T4 2023") +
  theme_void() +
  theme(legend.position = "top",
        legend.justification = "left",
        plot.title = element_text(margin = margin(b = 10)))

ggplot() +
  geom_sf(data = prov_esp, fill = "grey95", colour = "white") +
  geom_sf(data = paro_prov_sf, aes(size = total)) +
  scale_size(range = c(.5, 10)) + 
  labs(size = "Tasa", title = "Paro en el T4 2023") +
  theme_void() +
  theme(legend.position = "top",
        legend.justification = "left",
        plot.title = element_text(margin = margin(b = 10)))

ggplot() +
  geom_sf(data = prov_esp, fill = "grey90", colour = "white") +
  geom_sf(data = paro_prov_sf, aes(size = total, color = total)) +
  scale_size(range = c(.5, 10)) + # por área
  scale_color_binned_sequential(palette = "SunsetDark") +
  guides(color= guide_legend()) +
  labs(size = "Tasa", color = "Tasa", title = "Paro en el T4 2023") +
  theme_void() +
  theme(legend.position = "top",
        legend.justification = "left",
        plot.title = element_text(margin = margin(b = 10)))

Ejemplo: Mapa de accesibilidad media por carretera

En este caso importamos una geometría vectorial con las provincias de España, que contienen datos sobre la accesibilidad a los núcleos de población en 2015. Estos datos representan, en minutos, el promedio y la varianza de los tiempos de recorrido por carretera entre los núcleos de población (> 1500 habitantes por km2) de cada provincia.

prov_acces <- st_read("./data/accesibilidad.shp")
Reading layer `accesibilidad' from data source 
  `C:\Users\xeo19\Dropbox\RqueR\data\accesibilidad.shp' using driver `ESRI Shapefile'
Simple feature collection with 59 features and 9 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -18.15839 ymin: 27.63964 xmax: 4.31205 ymax: 43.75056
Geodetic CRS:  WGS 84
ggplot(prov_acces) +
        geom_sf(aes(fill = m))

Para cambiar el color de los límites provinciales a blanco definimos, dentro de la función geom_sf() el parámetro colour = "white".

# límites de países
world <- gisco_get_countries(resolution = "20")

ggplot(prov_acces) +
        geom_sf(aes(fill = var), 
                colour = "white")

# sin limitar
ggplot(prov_acces) +
        geom_sf(data = world) +
        geom_sf(aes(fill = var), 
                colour = "white")

p <- ggplot(prov_acces) +
        geom_sf(data = world, fill = "black", colour = "white") +
        geom_sf(aes(fill = var), 
                colour = "white") +
    coord_sf(xlim = c(-20, 5), 
             ylim = c(27, 45))
p

Aquí hacemos una pequeña paréntesis sobre dibujar límites internos de países. Aplicamos un con conjunto de funciones muy interesante a la ahora de simplificar geometría más.

# límites internos
wld_inner <- ms_innerlines(world)

plot(wld_inner)

# ahora con nuestros datos anteriores
p <- ggplot(prov_acces) +
        geom_sf(data = world, fill = "black", colour = NA) +
        geom_sf(data = wld_inner, 
                colour = "white",
                linewidth = .6) +
        geom_sf(aes(fill = var), 
                linewidth = .2,
                colour = "white") +
    coord_sf(xlim = c(-20, 5), 
             ylim = c(27, 45))

p

Y para cambiar la gama de colores de la variable usamos el paquete RColorBrewer, que contiene escalas de color predefinidas ampliamente utilizadas en cartografía (ver colorbrewer2.org). La función brewer.pal permite seleccionar el número de colores y la escala. Ponemos la leyenda en la parte inferior del gráfico con el parámetro legend.position = "bottom", consulta la ayuda de la función theme() y verás que tiene decenas de parámetros para modificar el aspecto del gráfico.

# paleta de colores
p <- p + scale_fill_distiller(palette = "RdYlGn") +
           labs(fill = "Tiempo (minutos)") +
              theme_bw() 
p

En ggplot también existe la posibilidad de modificar la leyenda con la función guides(). En ese caso tendríamos que usar los parametros de aes() (fill, size, colour, shape, etc.) para configurar la leyenda con datos continuos guide_colourbar() o con datos discretos guide_legend(). En la función guide_colourbar() se especifican varios elementos de estilo, desde la posición de las etiquetas hasta la anchura de la barra.

# cambiamos el estilo de la leyenda
p +  guides(fill = guide_colourbar()) +
          theme(legend.position = "bottom",
                legend.text.position = "bottom",
                legend.key.height = unit(.5, "lines"),
                legend.key.width = unit(3, "lines"),
                legend.title.position = "top")

# colorsteps
p <- p +  guides(fill = guide_colorsteps()) +
           theme(legend.position = "bottom",
                legend.text.position = "bottom",
                legend.key.height = unit(.5, "lines"),
                legend.key.width = unit(3, "lines"),
                legend.title.position = "top")

p

Al usar la clase sf no sólo podemos usar la función de geometría geom_sf(), sino también coord_sf(), que nos sirve para definir la proyección o limitar el área de visualización. Es decir, que es posible ajustar la proyección deseada en el proceso de la visualización sin tener que modificar los objetos espaciales originales.

# Cambiamos la proyección a LAEA
p + coord_sf(crs = 3035, 
             default_crs = 4326,
             xlim = c(-20, 5), 
             ylim = c(27, 45))
Coordinate system already present. Adding new coordinate system, which will
replace the existing one.

Por ejemplo podríamos proyectar, en el proceso de visualización, un mapa global sin la Antarctica a la proyección eurocéntrica LAEA (ETRS89 Lambert Azimuthal Equal-Area).

filter(world, NAME_ENGL != "Antarctica") |>
 ggplot() +
    geom_sf(fill = "black", 
            colour = "white",
            linewdith = .3) +
    coord_sf(crs = "epsg:3575") + 
    theme_light()
Warning in layer_sf(geom = GeomSf, data = data, mapping = mapping, stat = stat,
: Ignoring unknown parameters: `linewdith`

Cartografía más avanzada: el efecto hillshade

Es muy habitual ver mapas de relieve con efectos de sombras, también conocidos como ‘hillshade’ lo que genera una profundidad visual.

Paquetes

Paquete Descripción
tidyverse Conjunto de paquetes (visualización y manipulación de datos): ggplot2, dplyr, purrr,etc.
sf Simple Feature: importar, exportar y manipular datos vectoriales
elevatr Acceso a datos de elevación desde varias API
terra Importar, exportar y manipular raster (paquete sucesor de raster)
whitebox Una interfaz R para la biblioteca ‘WhiteboxTools’, que es una plataforma avanzada de análisis de datos geoespaciales
tidyterra Funciones auxilares para trabajar con {terra}
giscoR Límites administrativos del mundo
ggnewscale Extensión para ggplot2 de múltiples ‘scales’
ggblend Extensión para mezclar colores de gráficos ggplot
library(sf)
library(elevatr)
library(tidyverse)
library(terra)
library(ggnewscale)
library(tidyterra)
library(giscoR)
library(units)
library(ggblend)

Datos

Como área de interés usamos Suiza en este ejemplo. Con excepción de los límites de lagos, los datos necesarios los obtenemos a través de APIs usando diferentes paquetes. El paquete giscoR permite obtener los límites de países con diferentes resoluciones.

suiz <- gisco_get_countries(country = "Switzerland", resolution = "03")

plot(suiz)

Los límites de los lagos corresponden a una capa de modelos cartográficos digitales (DKM500) que ofrece swisstopo. El objetivo es quedar sólo con los grandes lagos, por tanto excluimos todos aquellos con menos de 50 km2 y también aquellos situados completamente en territorio italiano. Recuerda que con el paquete units podemos indicar unidades y así hacer cálculos.

# importamos los lagos
suiz_lakes <- st_read("./data/22_DKM500_GEWAESSER_PLY.shp")
Reading layer `22_DKM500_GEWAESSER_PLY' from data source 
  `C:\Users\xeo19\Dropbox\RqueR\data\22_DKM500_GEWAESSER_PLY.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 596 features and 14 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 2480000 ymin: 1062000 xmax: 2865000 ymax: 1302000
Projected CRS: CH1903+ / LV95
# filtramos a grandes lagos
suiz_lakes <- mutate(suiz_lakes, areakm = set_units(SHP_AREA, "m2") |> 
                                          set_units("km2")) |> 
                filter(areakm > set_units(50, "km2"),
                       !NAMN1 %in% c("Lago di Como / Lario",
                                     "Lago d'Iseo",
                                     "Lago di Garda"))
plot(suiz_lakes)
Warning: plotting the first 9 out of 15 attributes; use max.plot = 15 to plot
all

Modelo digital de terreno (MDT)

La función get_elev_raster() nos permite descargar un MDT de cualquier región del mundo a través de diferentes proveedores en formato de ráster. Por defecto usa AWS. Un argumento esencial es la resolución que depende de la latitud, la que se puede indicar como nivel de zoom (véase la ayuda). Por ejemplo, aquí usamos nivel 10 que a una latitud de 45º correspondería a aproximadamente 100m.

Después de obtener el MDT de Suiza debemos enmascarar los límites del país. La clase del objeto es RasterLayer del paquete raster, no obstante, el nuevo estándar es terra con la clase SpatRaster. Por eso lo convertimos y después aplicamos la máscara. Finalmente reproyectamos al sistema de coordenadas de Suiza obtenido de los datos vectoriales.

# obtenemos el mdt con 
mdt <- get_elev_raster(suiz, z = 7)
Mosaicing & Projecting
Note: Elevation units are in meters.
mdt # clase antigua de RasterLayer
class      : RasterLayer 
dimensions : 816, 1197, 976752  (nrow, ncol, ncell)
resolution : 0.004700506, 0.004700506  (x, y)
extent     : 5.625, 11.25151, 45.08689, 48.9225  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : file3e1411f34986.tif 
names      : file3e1411f34986 
plot(mdt)

# convertir a terra y enmascarar la zona de interés
mdt <- rast(mdt) |> 
         mask(vect(suiz)) 

# reproyectamos el mdt
mdt <- project(mdt, crs(suiz_lakes))

# reproyectamos los limite de lagos
suiz <- st_transform(suiz, st_crs(suiz_lakes))

Antes de calcular el efecto de sombra, creamos un simple mapa de relieve. En ggplot2 empleamos la geometría geom_raster() indicando la longitud, latitud y la variable para definir el color. Añadimos los límites de los lagos usando geom_sf() dado que se trata de un objeto sf. Aquí sólo indicamos el color de relleno con un azul claro. Con ayuda de scale_fill_hypso_tint_c() aplicamos una gama de colores correspondientes a un relieve, también llamado tintas hipsométricas, y definimos los cortes en la leyenda. En el resto de funciones hacemos ajustes de aspecto en la leyenda y en el estilo del gráfico.

# renombramos la capa de elevación
names(mdt) <- "alt"

# mapa
ggplot() +
  geom_spatraster(data = mdt,
              aes(fill = alt)) +
  geom_sf(data = suiz_lakes,
          fill = "#c6dbef", 
          colour = NA) +
  scale_fill_hypso_tint_c(breaks = c(180, 250, 500, 1000,
                                     1500,  2000, 2500,
                                     3000, 3500, 4000)) +
  guides(fill = guide_colorsteps()) +
  labs(fill = "m") +
  coord_sf() +
  theme_void() +
  theme(legend.position = "bottom",
        legend.key.height = unit(.5, "lines"),
        legend.key.width = unit(4, "lines"),
        legend.title.position = "right")
<SpatRaster> resampled to 501200 cells.

Internamente lo que hace geom_spatraster() es pasar el raster a un data.frame. geom_raster() sirve para una rejilla regular y geom_tile() si fuese irregular.

df <- as.data.frame(mdt, xy = T)

ggplot() +
  geom_raster(data = df,
              aes(x, y, fill = alt)) +
  theme_void()

Calcular el hillshade

Recordemos que el efecto hillshade es nada más que añadir una iluminación hipotética con respecto a una posición de una fuente de luz para así ganar profundidad. Las sombras dependen de dos variables, el acimut, el ángulo de la orientación sobra la superficie de una esfera, y la elevación, el ángulo de la altura de la fuente.

La información requerida para simular la iluminación es el modelo digital de terreno. La pendiente y la orientación podemos derivar del MDT usando la función terrain() del paquete terra. La unidad debe ser radianes. Una vez que tenemos todos los datos podemos hacer uso de la función shade() indicando el ángulo (elevación) y la dirección (acimut). El resultado es un ráster con valores entre 0 y 255, lo que nos indica sombras con bajos valores, siendo 0 negro y 255 blanco.

# estimamos la pendiente
sl <- terrain(mdt, "slope", unit = "radians")
plot(sl)

# estimamos la orientación
asp <- terrain(mdt, "aspect", unit = "radians")
plot(asp)

# calculamos el efecto hillshade con 45º de elevación 
hill_single <- shade(sl, asp, 
      angle = 45, 
      direction = 300,
      normalize= TRUE)

# resultado final hillshade 
plot(hill_single, col = grey(1:100/100))

Combinar el relieve y el efecto de sombra

El problema para añadir al mismo tiempo el relieve con su tinta hipsométrica y el efecto hillshade dentro de ggplot2 es que tenemos dos diferentes rellenos para cada capa. La solución consiste en usar la extensión ggnewscale que permite añadir múltiples scales. Primero añadimos con geom_raster() el hillshade, después definimos los tonos grises y antes de añadir la altitud incluimos la función new_scale_fill() para marcar otro relleno diferente. Para lograr el efecto es necesario dar un grado de transparencia a la capa del relieve, en este caso es del 70%. La elección de la dirección es importante, de ahí que debemos tener en cuenta siempre el lugar y el recorrido aparente del sol. (sunearthtools).

# nombre de la capa
names(hill_single)
[1] "hillshade"
# mapa 
ggplot() +
  geom_spatraster(data = hill_single,
              aes(fill = hillshade),
              show.legend = FALSE) +
  scale_fill_distiller(palette = "Greys", na.value = NA) +
  new_scale_fill() +
  geom_spatraster(data = mdt,
              aes(fill = alt),
              alpha = .7) +
  scale_fill_hypso_tint_c(breaks = c(180, 250, 500, 1000,
                                     1500,  2000, 2500,
                                     3000, 3500, 4000)) +
  geom_sf(data = suiz_lakes,
          fill = "#c6dbef", colour = NA) +
  guides(fill = guide_colorsteps()) +
  labs(fill = "m") +
  coord_sf() +
  theme_void() +
    theme(legend.position = "bottom",
        legend.key.height = unit(.5, "lines"),
        legend.key.width = unit(4, "lines"),
        legend.title.position = "right")
<SpatRaster> resampled to 501200 cells.
<SpatRaster> resampled to 501200 cells.

Sombras multidireccionales

Lo que hemos visto es un efecto unidireccional, aunque es lo más habitual, podemos crear un efecto más suave e incluso más realista combinando varias direcciones.

Simplemente mapeamos sobre un vector de varias direcciones al que se aplica la función shade() con una elevación fija. Después convertimos la lista de ráster en un objeto multidimensional de varias capas para reducirlas sumando todas las capas.

# pasamos varias direcciones a shade()
hillmulti <- map(c(270, 15, 60, 330), function(dir){ 
                    shade(sl, asp, 
                          angle = 45, 
                          direction = dir,
                          normalize= TRUE)}
  )

# creamos un raster multidimensional y lo reducimos sumando
hillmulti <- rast(hillmulti) |> sum()

# multidireccional
plot(hillmulti, col = grey(1:100/100))

# unidireccional
plot(hill_single, col = grey(1:100/100))

Hacemos lo mismo como antes para visualizar el relieve con sombras multidireccionales.

# convertimos el hillshade a xyz
names(hillmulti) <- "hillshade"

# mapa
ggplot() +
  geom_spatraster(data = hillmulti,
              aes(fill = hillshade),
              show.legend = FALSE) +
  scale_fill_distiller(palette = "Greys", na.value = NA) +
  new_scale_fill() +
  geom_spatraster(data = mdt,
              aes(fill = alt),
              alpha = .7) +
  scale_fill_hypso_tint_c(breaks = c(180, 250, 500, 1000,
                                     1500,  2000, 2500,
                                     3000, 3500, 4000)) +
  geom_sf(data = suiz_lakes,
          fill = "#c6dbef", colour = NA) +
  guides(fill = guide_colorsteps()) +
  labs(fill = "m") +
  coord_sf() +
  theme_void() +
        theme(legend.position = "bottom",
        legend.key.height = unit(.5, "lines"),
        legend.key.width = unit(4, "lines"),
        legend.title.position = "right")
<SpatRaster> resampled to 501200 cells.
<SpatRaster> resampled to 501200 cells.

La técnica de mezcla de colores es muy útil para obtener resultados notables en el efecto de sombreado. Desde hace poco el paquete ggblend ofrece esta posibilidad. Con el objetivo de combinar varias capas, es necesario insertar los objetos geom_raster() y los scale_fill_*() en una lista seperados por coma. Después le sigue el pipe con la función blend("tipo_de_mezcal") al que le sumamos los otros objetos de ggplot2. En este caso aplicamos la multiplicación como forma de mezcla.

# mapa
m <- ggplot() +
  list(
  geom_spatraster(data = hillmulti,
              aes(fill = hillshade),
              show.legend = FALSE),
  scale_fill_distiller(palette = "Greys", na.value = NA),
  new_scale_fill(),
  geom_spatraster(data = mdt,
              aes(fill = alt),
              alpha = .7),
  scale_fill_hypso_tint_c(breaks = c(180, 250, 500, 1000,
                                     1500,  2000, 2500,
                                     3000, 3500, 4000))
  ) |> blend("multiply") +
  geom_sf(data = suiz_lakes,
          fill = "#c6dbef", colour = NA) +
  guides(fill = guide_colorsteps()) +
  labs(fill = "m") +
  coord_sf() +
  theme_void() +
      theme_void() +
        theme(legend.position = "bottom",
        legend.key.height = unit(.5, "lines"),
        legend.key.width = unit(4, "lines"),
        legend.title.position = "right")

ggsave("mdt_hillshade_blend.png", m, 
       width = 10, 
       height = 8, 
       unit = "in",
       device = png, 
       type = "cairo",
       bg = "white")